Load libraries

library(tidyverse)
library(scran)
library(scater)
library(gprofiler2)
library(gridExtra)
library(ComplexHeatmap)
library(circlize)
library(edgeR)
library(readxl)
library(ggrepel)
library(tftargets)
library(network)
library(ggnet)
library(parallel)
library(BiocParallel)
library(ggsci)
library(airr)
library(numbat)
library(scales)

ncores <- 10
mcparam <- MulticoreParam(workers=ncores)
register(mcparam)

Define variables

opt <- list()
opt$pathinfo <- "misc/pathways_2.xlsx"
opt$scesc <- "data/submission/scRNA_complete.RDS"
opt$scescpat <- "data/submission/scRNA_per_pat_complete.RDS"
opt$scebulk <- "data/submission/RNA_complete.RDS"
opt$plot <- "plots/SFig4_5/"


## Set theme
lgd <-  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), 
        text = element_text(size = 15),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA), 
        legend.background = element_rect(fill='transparent',colour = NA), 
        strip.background = element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

Load data

## SCE
sce <- readRDS(opt$scebulk)
sce
## class: SingleCellExperiment 
## dim: 36601 117 
## metadata(0):
## assays(2): counts logcounts
## rownames(36601): ENSG00000000003 ENSG00000000005 ... ENSG00000288459
##   ENSG00000288460
## rowData names(3): ID Symbol Type
## colnames(117): H358_BafilomycinA1 H358_Birinapant ... H524_Selinexor
##   H524_Venetoclax
## colData names(5): uniqueBarcode patientID Diagnosis_simple Treatment
##   label
## reducedDimNames(2): PCA UMAP
## mainExpName: NULL
## altExpNames(0):
## scRNA-Seq data
sce.sc <- readRDS(opt$scesc)

## single-cell cell types
types <- sce.sc[["types"]]

## Load tcr data
immTab <- sce.sc[["immTab"]]

sce.sc <- sce.sc[["SCE"]]

## sce per pat
sce.pat <- readRDS(opt$scescpat)

## Load pathway annotation
path.info <- read_excel(opt$pathinfo)

Load transcription factors

ntop <- 1000
maxtargets <- 5000

## Get TF target numbers
motif.df <- lapply(names(Marbach2016), function(x) {
  tf.list <- Marbach2016[[x]]
  #tf.list <- TRRUST[[x]]
  
  data.frame(motif_name = x, n_targets = length(tf.list))
  
}) %>% bind_rows()

## Get TF targets
target.df <- lapply(names(Marbach2016), function(x) {
  tf.list <- Marbach2016[[x]]
  
  data.frame(motif_name = x, targets = tf.list)
  
}) %>% bind_rows()

Analysis

Perform differential gene expression testing

testTreat <- unique(sce$Treatment)

de.res <- mclapply("DMSO", mc.cores = ncores, 
                 function(z) {
                   
  message(paste0("Fitting model for ", z))
                   
  compDrug <- testTreat[testTreat != z] # instead of subtypes we use other treatments here
                   
  mclapply(compDrug, mc.cores = ncores, 
         function(x) {
    message(paste0("Drug ", x))
    summed <- sce#[hvgenes, ]
    colData(summed)$Treatment <- factor(colData(summed)$Treatment, levels = c(z, compDrug))
    y <- DGEList(counts(summed), samples = colData(summed))
    design <- model.matrix(~0 + factor(patientID) + Treatment, y$samples) # set coefficient (name or matrix column) - fig. 7 f1000
    keep <- filterByExpr(y, design = design)
    y <- y[keep,]
    summary(keep)
    message(paste0("Keeping ", sum(keep), " genes"))
    y <- calcNormFactors(y)
    y <- estimateDisp(y, design) # gene-wise dispersion (variance on top of poisson distribution)
    plotBCV(y)
    fit <- glmQLFit(y, design, robust = TRUE) # QL way of fitting the model
    plotQLDisp(fit)
    res <- glmQLFTest(fit, coef = paste0("Treatment", x))
    summary(decideTests(res))
    resTab <- res$table
    resTab$gene <- rownames(resTab)
    resTab <- resTab %>% mutate(p.adj = p.adjust(PValue, method = "BH")) #%>%
      #mutate(p.adj = ifelse(p.adj < 1e-12, 1e-12, p.adj))
    resTab$Treatment <- x
    resTab$Control <- z
    resTab %>% arrange(desc(logFC))
  }) %>% bind_rows()
}) %>% bind_rows()
## Fitting model for DMSO
de.res <- de.res %>%
  left_join(select(data.frame(rowData(sce)), ID, Symbol), c("gene" = "ID"))

#write.csv(de.res, opt$deres)

sigTab <- de.res %>%
  group_by(Treatment) %>%
  filter(p.adj < 0.1) %>%
  filter(logFC > 0.25 | logFC < -0.25) %>%
  mutate(dir = ifelse(logFC > 0, "Up", "Down")) %>%
  dplyr::count(dir) %>%
  dplyr::rename(Direction = dir) %>%
  ungroup()

drugOrder <- sigTab %>%
  group_by(Treatment) %>%
  summarise(allDEG = sum(n)) %>%
  arrange(desc(allDEG)) %>%
  pull(Treatment)

Construct a network of differentially expressed genes

# Assign nodes (treatments) and edges (5% FDR DEGs)
plotMat <- de.res %>% filter(p.adj < 0.1) %>%
  filter(logFC > 0.25 | logFC < -0.25)
allTreat <- unique(plotMat$Treatment)
target <- data.frame(Treatment = allTreat,
                     target = c("Autophagyi", "IAPi", "HDACi", "mTORi", "ITKi", "TLR8a", "MDM2i", 
                                "IMID", "JAKi", "XPOi", "BCL2i"))

opt$sub <- TRUE

if(opt$sub == TRUE) {
  plotMat <- de.res %>% filter(p.adj < 0.1, !Treatment %in% c("Pomalidomide", "Dacinostat")
                             ) %>%
  filter(logFC > 0.25 | logFC < -0.25) 
} else {
  plotMat <- de.res %>% filter(p.adj < 0.1) %>%
    filter(logFC > 0.25 | logFC < -0.25)
}


allTreat <- unique(plotMat$Treatment)

nodes <- data.frame(Treatment = allTreat)
edges <- plotMat %>% select(Treatment, Symbol) %>%
  mutate(Treatment = as.character(Treatment))

# Create network object
set.seed(100)
net <- network(edges, directed = TRUE)
netTab <- ggnet2(net)$data |> 
  mutate(Type = ifelse(label %in% allTreat, "Treatment", "Symbol"),
         Treatment = ifelse(Type == "Treatment", label, NA)) |> 
  select(-alpha,-color,-shape,-size) |> 
  mutate(x = as.numeric(x), y = as.numeric(y))
  
# Define edges based on direction and padj of DEGs
from <- edges %>%
  left_join(netTab, by = c("Treatment" = "label")) %>% dplyr::select(Treatment, x, y) %>%
  dplyr::rename(x_start = x, y_start = y) |> distinct()

to <- edges %>%
  left_join(netTab, by = c("Symbol" = "label")) %>% dplyr::select(Symbol,x,y) %>%
  dplyr::rename(x_end = x, y_end = y) |> distinct()

edges <- edges %>% left_join(from) %>% left_join(to) %>%
  left_join(plotMat %>% select(Symbol,Treatment,logFC,p.adj)) %>%
  distinct() %>%
  mutate(Direction = ifelse(logFC < 0,"Down","Up"),
         alpha = -log10(p.adj)
         #alpha = rescale(-log10(p.adj#squish(p.adj, c(1.812407e-209,1))))
         )

# Create auxiliary table for annotation

n_DEG <- plotMat %>% #filter(!Treatment %in% c("Pomalidomide", "BafilomycinA1", "Dacinostat")) %>%
  group_by(Treatment) %>%
  filter(p.adj < 0.1) %>%
  filter(logFC > 0.25 | logFC < -0.25) %>% 
  dplyr::count() %>% ungroup()

netTabAnno <- netTab |> 
  left_join(target) |> 
  left_join(n_DEG) |> 
  mutate(Treatment = factor(Treatment, levels = allTreat))

p <- ggplot()  +
  geom_curve(data = edges,
               aes(x = x_start, xend = x_end, y = y_start, yend = y_end, col = Direction, alpha = p.adj),
               linewidth = 0.1, curvature = 0.05) +
  geom_point(data = netTabAnno[netTabAnno$Type == "Symbol",],
             aes(x = x, y = y),
             size = 0.5, shape = 21, color = "grey60", fill = "grey95", stroke = 0.75, alpha = 0.6) +
  geom_point(data = netTabAnno[netTabAnno$Type == "Treatment",],
             aes(x = x, y = y, fill = Treatment, size = n), shape = 21, stroke = 0.75) +
  # scale_fill_manual(values = col_TX) +
  scale_size_continuous(name = "Number of DEGs",
                        breaks = seq(0, max(n_DEG$n), by = 1000), limits = c(0, max(n_DEG$n)), range = c(0.5,7.5)) +
  scale_color_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
  scale_alpha_continuous() +
  geom_label_repel(data = netTabAnno |> filter(Type == "Treatment"),
                  aes(x = x, y = y, label = Treatment, fill = Treatment, point.size = rescale(n, to = c(0.5,7.5))),
                  max.overlaps = Inf, vjust = "bottom", nudge_y = -0.001, fontface = "bold", alpha = 0.7, #direction = "y", 
                  show.legend = FALSE, size = 5)  +
  guides(color = guide_legend(ncol = 2, override.aes = list(linewidth=2.5)),
         size = guide_legend(ncol = 2),
         alpha = guide_legend(override.aes = list(linewidth=2.5)),
         fill = guide_legend(override.aes = list(size=3.5))) +
  xlab("") + ylab("") + ggtitle("Network of Differentially Expressed Genes") +
  lgd +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        panel.border = element_rect(linewidth = 1, colour = NA, fill = "transparent")) +
  scale_fill_manual(values = c("Venetoclax" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000", 
                                "Ibrutinib" = "#283593", "Everolimus" = "#43A047", "Selinexor" = "#FFD45F",
                               "BafilomycinA1" = "#F8BBD0", "Motolimod" = "#AB47BC", "Dacinostat" = "#F44336", 
                                "Dacinostat" = "#64B5F6", "Pomalidomide" = "#80DEEA", "Ruxolitinib" = "#4DB6AC"))

p

if(opt$sub == TRUE) {
 ggsave(plot = p, filename = paste0(opt$plot, "deg_subnetwork.png"),
       height = 20, width = 25, units = "cm")
} else {
  ggsave(plot = p, filename = paste0(opt$plot, "deg_network.png"),
       height = 20, width = 25, units = "cm")
}

Volcano plots

pList <- lapply(c("Motolimod"), function(x) {
  
  if(x %in% c("Birinapant", "Motolimod", "Everolimus", "Ibrutinib", "BaflimoycinA1")) {
    geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
  } else if(x %in% c("Nutlin-3a", "Selinexor")) {
    geneSelec <- path.info %>% filter(pathway %in% c("TP53up")) %>% select(gene) %>% unlist() %>% as.vector()
    geneSelec <- c("TP53", geneSelec)

  } else if(x %in% c("Everolimus")) {
    geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
  }
  
  volTab <- de.res %>% filter(Treatment == x)

  ## Define limits  
  minLim <- -1 * max(abs(de.res$logFC))
  maxLim <- max(abs(de.res$logFC))
  thresh <- 0.25

  volTab$diffCol <- NA
  volTab[volTab$logFC > 0 & volTab$p.adj < 0.1 & volTab$logFC > thresh, "diffCol"] <- "sigUp"
  volTab[volTab$logFC < 0 & volTab$p.adj < 0.1 & volTab$logFC < -thresh, "diffCol"] <- "sigDown"
  volTab[volTab$p.adj >= 0.1, "diffCol"] <- "ns"
  volTab[volTab$logFC < thresh & volTab$logFC > -(thresh), "diffCol"] <- "ns"
  
  p <- volTab %>%
    mutate(logFC = ifelse(logFC > maxLim, maxLim, logFC)) %>%
    mutate(logFC = ifelse(logFC < minLim, minLim, logFC)) %>%
    #filter(gene %in% geneSelec) %>%
    filter(Treatment == x) %>%
    ggplot(aes(x = logFC, y = -log10(p.adj), fill = diffCol)) +
    geom_point(shape = 21, size = 2) + xlab("Log2FC") + ylab("-Log10(p.adj)") + ggtitle(paste0(x)) +
    #xlim(minLim, maxLim) +
    geom_hline(yintercept = -log10(0.1), linetype = "dashed") +
    geom_vline(xintercept = -thresh, linetype = "dashed") + geom_vline(xintercept = thresh, linetype = "dashed") +
    geom_label_repel(data = filter(volTab, Symbol %in% c(geneSelec), Treatment %in% x, p.adj < 0.1, 
                                 logFC > thresh | logFC < -(thresh)), aes(label = Symbol), fill = "white", 
                    nudge_y = 0.25, size = 4, max.overlaps = getOption("ggrepel.max.overlaps", default = 30), 
                     segment.linetype = 2, force = 20, alpha = 0.8) +
    scale_fill_manual(values = c("sigUp" = "#F44336", "sigDown" = "#1976D2", "ns" = "grey")) +
    theme_classic() +
    theme(plot.title = element_text(hjust = 0.5), 
          legend.position = "none", 
          text = element_text(size = 17.5), 
          axis.text = element_text(size = 15),
          panel.background = element_rect(fill = "transparent",colour = NA),
          plot.background = element_rect(fill = "transparent",colour = NA),
          legend.background = element_rect(fill='transparent'),
          strip.background = element_blank())
  
   ggsave(plot = p, filename = paste0(opt$plot, "vol_", x, ".png"),
          height = 15, width = 15, unit = "cm")
  p

})
pList 
## [[1]]

#### Plot bafilomycin A1 vs Motolimod
p <- de.res %>% filter(Treatment %in% c("BafilomycinA1", "Motolimod")) %>%
  select(gene, Symbol, logFC, Treatment) %>%
  pivot_wider(names_from = "Treatment", values_from = logFC) %>%
  ggplot(aes(x = BafilomycinA1, y = Motolimod)) +
  xlab("Bafilomycin A1 (log2FC)") + ylab("Motolimod (log2FC)") +
  ggtitle("Motolimod vs Bafilomycin A1") +
  geom_smooth(method = "lm", colour = "black", fill = "grey80") +
  geom_point() +
    theme_classic() +
    theme(plot.title = element_text(hjust = 0.5), 
          legend.position = "none", 
          text = element_text(size = 17.5), 
          axis.text = element_text(size = 15),
          panel.background = element_rect(fill = "transparent",colour = NA),
          plot.background = element_rect(fill = "transparent",colour = NA),
          legend.background = element_rect(fill='transparent'),
          strip.background = element_blank())

p
## `geom_smooth()` using formula = 'y ~ x'

ggsave(plot = p, filename = paste0(opt$plot, "moto_vs_bafilo.png"),
       height = 15, width = 15, units = "cm")
## `geom_smooth()` using formula = 'y ~ x'

TP53 dot plot

geneSelec <- path.info %>% filter(pathway %in% c("TP53up")) %>% select(gene) %>% unlist() %>% as.vector()

geneSub <- lapply(unique(de.res$Treatment), function(x) {
  de.res %>% 
    filter(Treatment == x) %>%
    filter(p.adj < 0.1) %>%
    filter(logFC > 0.25 | logFC < -0.25) %>%
    pull(Symbol)
  
}) %>% unlist() %>% unique()

geneSub <- intersect(geneSelec, geneSub)

plotMat.wide <- de.res %>% 
  mutate(logFC = ifelse(p.adj >= 0.1, 0, logFC)) %>%
  dplyr::select(Symbol, Treatment, logFC) %>% 
  filter(Symbol %in% geneSub) %>%
  dplyr::rename(log2FC = logFC) %>%
  pivot_wider(names_from = "Treatment", values_from = "log2FC") %>%
  column_to_rownames("Symbol")

geneOrder <- rownames(plotMat.wide)[hclust(dist(plotMat.wide), method = "ward.D2")$order]


p <- de.res %>% filter(#Treatment %in% c("Birinapant", "Motolimod"),
  Symbol %in% geneSub) %>%
  #filter(Symbol == "MDM2") %>%
  filter(p.adj < 0.1) %>%
  mutate(pSign = -log10(p.adj) * sign(logFC),
         p.adj = ifelse(-log10(p.adj) > 5, 5, -log10(p.adj)), 
         Symbol = factor(Symbol, levels = geneOrder)) %>%
  dplyr::rename(log2FC = logFC) %>%
  ggplot(aes(y = Symbol, x = Treatment)) +
  geom_point(aes(size = p.adj, fill = log2FC), shape = 21) +
  xlab("") + ylab("") + ggtitle("Effect on TP53 Target Genes") +
  lgd +
  scale_fill_gradient2(high = "#F44336", low = "#1976D2", name = "log2FC",
                       limits=c(-2, 2), oob=squish) +
  scale_size_continuous(name = "-Log10(p.adj)",
                        breaks = seq(0, 5, by = 2.5),
                        limits = c(0, 5)
                        ) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
        text = element_text(size = 12.5),
        panel.border = element_rect(linewidth = 1))
p

ggsave(plot = p, filename = paste0(opt$plot, "tp53_long.png"),
       height = 40, width = 12.5, units = "cm")

GSEA

go.df.merge <- mclapply(c("Birinapant", "Motolimod", "Selinexor", "BafilomycinA1", "Dacinostat", "Pomalidomide", "Ibrutinib", 
                        "Everolimus", "Venetoclax", "Ruxolitinib", "Nutlin-3a"), mc.cores = ncores, 
                      function(compTreat) {
    message(paste("Running gprofiler for", compTreat, sep = " "))
    sigUp <- de.res %>% filter(logFC > 0.25, Treatment == compTreat, 
                               #!gene %in% c(drugOther), 
                               p.adj <= 0.1) %>% arrange(desc(logFC))
    #sigUp <- de.res %>% filter(celltype == compType, logFC > 0.25, Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
    
    sigUp
    if(nrow(sigUp) > 1) {
      gostresUp <- gost(query = sigUp$gene,
                        organism = "hsapiens", ordered_query = TRUE,
                        multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
                        measure_underrepresentation = FALSE, evcodes = FALSE,
                        user_threshold = 0.1, correction_method = "fdr",
                        domain_scope = "annotated", custom_bg = NULL,
                        numeric_ns = "", sources = NULL, as_short_link = FALSE)
      gostresUp$result

      go.df.up <- gostresUp$result %>% arrange(p_value)
      go.df.up$rank <- 1:nrow(go.df.up)
      go.df.up$type <- "Up"
    } else {
      go.df.up <- data.frame(p_value = NA, rank = NA, type = "Up", source = NA, term_name = NA, 
                               intersection_size = 10)
    }

    sigDown <- de.res %>% filter(logFC < -(0.25), Treatment == compTreat,
                                 #!gene %in% c(drugOther),
                                  p.adj <= 0.1) %>% arrange(logFC)
    #sigDown <- de.res %>% filter(celltype == compType, logFC < -(0.25), Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
    if(nrow(sigDown) > 1) {
    gostresDown <- gost(query = sigDown$gene,
                        organism = "hsapiens", ordered_query = TRUE,
                        multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
                        measure_underrepresentation = FALSE, evcodes = FALSE,
                        user_threshold = 0.1, correction_method = "fdr",
                        domain_scope = "annotated", custom_bg = NULL,
                        numeric_ns = "", sources = NULL, as_short_link = FALSE)

    gostresDown$result

    go.df.down <- gostresDown$result %>% arrange(p_value)
    go.df.down$rank <- 1:nrow(go.df.down) * -1
    go.df.down$type <- "Down"
    } else {
      go.df.down <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA, 
                               intersection_size = 10)
    }

    go.df.merge <- rbind(go.df.up, go.df.down)
    go.df.merge$Treatment <- compTreat

    go.df.merge <- go.df.merge %>% filter(!intersection_size < 10) %>% 
      filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
      mutate(p.adj = p.adjust(p_value, method = "BH"))
    
    if(nrow(go.df.merge) < 1) {
      go.df.merge <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA, 
                               intersection_size = 0)
    } else {
      go.df.merge
    }
}) %>% bind_rows() 

## Exclude all terms that are up and downregulated by a single-drug
dblTerm <- go.df.merge %>% 
  filter(!intersection_size < 10) %>%
  filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
  mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
  #filter(term_name %in% unique(termTab$term_name)) %>%
  dplyr::count(term_name, Treatment) %>%
  filter(n > 1) %>%
  pull(term_name) %>% unique()

plotMat <- go.df.merge %>% 
  mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
  filter(!source %in% c("CORUM", "TF", "MIRNA"), !intersection_size < 10) %>%
  filter(!term_name %in% c(dblTerm)) %>%
  #filter(term_name %in% gseaSub) %>%
  select(Treatment, term_name, type, p.adj) %>%
  mutate(p.adj = -log10(p.adj), 
         p.adj = ifelse(type == "Down", p.adj * -1, p.adj)) %>%
  select(-type) %>%
  pivot_wider(names_from = "term_name", values_from = "p.adj") %>%
  column_to_rownames("Treatment")

plotMat[is.na(plotMat)] <- 0

pHeat <- Heatmap(plotMat, show_column_names = FALSE, 
        clustering_method_rows = "ward.D2", 
        clustering_method_columns = "ward.D2", border = TRUE, 
        column_title = paste0("Enriched Gene Sets Among DEGs (FDR < 10%)"), 
        column_title_gp = gpar(fontsize = 17.5),
        heatmap_legend_param = list(title_gp = gpar(fontsize = 12.5)),
        col=colorRamp2(c(-5, 0, 5), colors = c("#1976D2", "white", "#F44336")), 
        name = "-Log10(p.adj)\nwith Direction")
## Warning: The input is a data frame-like object, convert it to a matrix.
pHeat

png(file=paste0(opt$plot, "GSEA_heatmap.png"),
    height = 10, width = 25, res = 600, units = "cm", bg = "transparent")
draw(pHeat, background = "transparent")
dev.off()
## png 
##   2
## Plot top enriched pathwaysgo.df.merge
nTop <- 20

lapply(c("Selinexor"), function(y) {
           
  gseaTab <- go.df.merge %>% filter(Treatment == y) %>% filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
        mutate(term_name = paste0(source, "_", term_name)) %>% unique()
  
  message(y)
     lapply(unique(gseaTab$type), function(z) {
      message(z)
      
      gseaTab <- gseaTab %>% filter(type == z)
      
      gseaTab <- gseaTab[1:nTop, ] # get only top enrichments
      
      termOrder <- gseaTab %>% arrange(desc(p_value)) %>% pull(term_name)
      
      gseaTab <- gseaTab %>%
        mutate(term_name = factor(term_name, levels = termOrder), 
               p.adj = p.adjust(p_value, method = "BH"))
      
    if(-log10(min(gseaTab$p.adj, na.rm = TRUE)) > 10) {
      x.max <- -log10(min(gseaTab$p.adj, na.rm = TRUE))
    } else {
      x.max <- 10
    }
    if(y %in% c("Birinapant", "Selinexor", "Pomalidomide")) {
      x.max <- 3
    }
    if(y %in% c("Pomalidomide")) {
      x.max <- 2
    }

      
      p <- gseaTab %>%
        ggplot(aes(x = -log10(p.adj), y = term_name)) +
        geom_bar(stat = "identity", colour = "black", alpha = 0.5, aes(fill = type)) +
        geom_vline(xintercept = -log10(0.1), alpha = 0.8, linetype = "dashed") +
        geom_text(aes(label = term_name, y = term_name, x = 0.25), hjust = 0, size = 4.25) +
        xlab("-Log10(p.adj)") + ylab("") + #ggtitle(paste0("Top enriched gene-set: ", z)) +
        ggtitle(paste0(y, " - ", z)) +
        xlim(0, x.max) +
        lgd + theme(axis.text.y = element_blank(), 
                    axis.ticks.y = element_blank(), 
                    text = element_text(size = 22.5),
                    panel.border = element_rect(colour = "black", fill=NA, size=1.5), 
                    strip.background = element_blank(),
                    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                    panel.background = element_rect(fill = "transparent",colour = NA),
                    plot.background = element_rect(fill = "transparent",colour = NA),
                    legend.background = element_rect(fill='transparent'), 
                    legend.position = "none") +
        scale_fill_manual(values = c("Down" = "#1976D2", "Up" = "#F44336"))
        #scale_fill_npg() +
  #guides(fill=guide_legend(title="Source"))
      ggsave(plot = p,
             file = paste0(opt$plot, "GSEA/", y, "_", z, "_top_", nTop, ".png"),
             height = 22.5, width = 19, units = "cm")
      p
    })
})
## [[1]]
## [[1]][[1]]

Hypergeometric test for transcription factors

treatSub <- c("Pomalidomide")

## Run hypergeometric test
hyper.df <- lapply(treatSub, function(x) {
  ## Limit the number of targets
  motif.df <- motif.df %>% filter(n_targets < maxtargets)
  message("Keeping ", nrow(motif.df), " TFs with less than ", maxtargets)

  ## Get the transcription factor targets
  lapply(c("Up", "Down"), function(y) {
    df.hub <- de.res %>% filter(p.adj < 0.05, Treatment == x)
    
    if(y == "Up") {df.hub <- df.hub %>% filter(logFC > 0.25)}
    if(y == "Down") {df.hub <- df.hub %>% filter(logFC < -0.25)}
    
    mclapply(unique(motif.df$motif_name), mc.cores = ncores, 
             function(z) {
      message(z)
      tf.targets <- target.df %>% filter(motif_name == z) %>% pull(targets)

      size.module <- nrow(df.hub)
      size.genome <- length(unique(de.res$Symbol))
      size.genome <- nrow(sce)
      n.targets <- length(tf.targets)
      n.over <- intersect(df.hub$Symbol, tf.targets) %>% length()

      p.val <- phyper(q = n.over - 1, # number of white balls drawn
                      m = n.targets, # number white balls in urn
                      n = size.genome - n.targets, # number black balls in urn
                      k = size.module, # number of draws
                      lower.tail= FALSE)

      df <- data.frame(Treatment = x, dir = y, tf = z, p.val = p.val, n.over = n.over,
                       n.targets = n.targets, size.module = size.module,
                       size.genome = size.genome)
      #df %>% filter(n.over != 0)
    }) %>% bind_rows()
  }) %>% bind_rows() %>% mutate(p.adj = p.adjust(p.val, method = "BH"))
})# %>% bind_rows()

names(hyper.df) <- treatSub

lapply(c("Pomalidomide"), function(x) {
  lapply(c("Up", "Down"), function(y) {
    message(x)
    hyper.df <- hyper.df %>% bind_rows() %>%
      filter(Treatment == x, dir == y) %>%
      arrange(p.adj) %>%
      filter(p.adj < 0.1) %>%
      mutate(p.adj = -log10(p.adj))

    hyper.df$rank <- 1:nrow(hyper.df)

    ## Get the top 5 enriched TPs
    toptf <- hyper.df %>% filter(rank <= 10) %>%
      pull(tf)

    if(x == "Pomalidomide") {
      toptf <- c(toptf, "IKZF1", "IKFZ3", "BRD7", "ARID2", "IRF4", "MYC")
    }

    if(x %in% c("Nutlin-3a", "Selinexor")) {
      toptf <- c(toptf, "TP53")
    }
    if(x %in% c("Birinapant", "Motolimod", "BafilomycinA1", "Ibrutinib")) {
      toptf <- c(toptf, "NFKB", "NFKB1", "NFKB2", "REL", "RELA", "RELB")
    }

    if(y == "Up") {col.point <- "#F44336"}
    if(y == "Down") {col.point <- "#1976D2"}

    p <- hyper.df %>%
      ggplot(aes(x = rank, y = p.adj)) +
      geom_point() +
      geom_point(data = hyper.df[hyper.df$tf %in% toptf, ], fill = col.point, size = 2.5, shape = 21) +
      xlab("Rank") + ylab("-Log10(p.adj)") + ggtitle(paste0(x, " - TF Enrichment (", y, ")")) +
      geom_label_repel(aes(label = tf), data = hyper.df[hyper.df$tf %in% toptf, ],
                       max.overlaps = getOption("ggrepel.max.overlaps", default = 30),
                       segment.linetype = 2, force = 20, alpha = 0.8) +
      lgd +
      theme_classic() +
      theme(plot.title = element_text(hjust = 0.5),
            legend.position = "none",
            text = element_text(size = 17.5),
            axis.text = element_text(size = 15),
            panel.background = element_rect(fill = "transparent",colour = NA),
            plot.background = element_rect(fill = "transparent",colour = NA),
            legend.background = element_rect(fill='transparent'),
            strip.background = element_blank())

    ggsave(plot = p, filename = paste0(opt$plot, x, "_", y, "_tf_enrichment.png"),
           height = 12.5, width = 16, units = "cm")

    p
  })
})
## [[1]]
## [[1]][[1]]

## 
## [[1]][[2]]

Overview scRNA-Seq data

## Randomly sample 
set.seed(100)
sce.sc <- sce.sc[, sample(1:ncol(sce.sc), ncol(sce.sc))]

## Re-run UMAP
set.seed(100)
sce.sc <- runUMAP(sce.sc, n_neighbors = 15, min_dist = 0.75, name = "UMAP", BPPARAM = mcparam)

## Make plots
plotTab <- ggcells(sce.sc)[[1]]

## Pairwise each treatment
lapply(c("Birinapant", "Nutlin-3a", "Ibrutinib", "Everolimus", "Selumetinib"), function(x) {
  p <- plotTab %>%
    filter(Treatment %in% c("DMSO", x)) %>%
    ggplot(aes(x = UMAP.1, y = UMAP.2, fill = Treatment)) +
    geom_point(shape = 21, size = 2) +
    lgd + theme_void() + theme(legend.position = "none") + 
    scale_fill_manual(values = c("DMSO" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000", 
                                 "Ibrutinib" = "#283593", "Everolimus" = "#43A047", 
                                 "Selumetinib" = "#F8BBD0"))
  
    ggsave(plot = p, filename = paste0(opt$plot, x, "_scrna_umap.png"), 
         height = 15, width = 15, units = "cm")
    p
})
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Show donut plots

plotTab.vdj <- plotTab %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(celltype %in% c("malignant B-cell", "healthy B-cell")) %>%
  mutate(celltype = ifelse(celltype == "malignant B-cell", patientID, celltype), 
         celltype = ifelse(celltype == "healthy B-cell", "Healthy B-Cells", celltype))

plotTab.vdj[plotTab.vdj$celltype == "H525", "celltype"] <- "HCL1"
plotTab.vdj[plotTab.vdj$celltype == "H526", "celltype"] <- "HCL2"
plotTab.vdj[plotTab.vdj$celltype == "H432", "celltype"] <- "MCL1"
plotTab.vdj[plotTab.vdj$celltype == "H452", "celltype"] <- "MCL2"
plotTab.vdj[plotTab.vdj$celltype == "H496", "celltype"] <- "MCL3"

## Merge VDJ data with plotting info
plotTab.vdj <- plotTab.vdj %>% left_join(select(immTab, cell_id, clone_id, c_call#, mu_freq
                                                ), 
                                     by = c("uniqueBarcode" = "cell_id"))

plotTab.vdj %>% filter(clone_id %in% c("60_5")) %>% select(celltype)
## [1] celltype
## <0 rows> (or 0-length row.names)
## Calculate the number of clones per cell type
countTab <- plotTab.vdj %>%
  group_by(celltype) %>%
  select(celltype, clone_id) %>% unique() %>%
  dplyr::count(celltype) %>% mutate(n_clones = n) %>%
  select(-n)
countTab
## # A tibble: 6 × 2
## # Groups:   celltype [6]
##   celltype        n_clones
##   <chr>              <int>
## 1 HCL1                   2
## 2 HCL2                   2
## 3 Healthy B-Cells      176
## 4 MCL1                   2
## 5 MCL2                   2
## 6 MCL3                   2
## Calculate the size of each clone
countSize <- plotTab.vdj %>%
  filter(!is.na(clone_id)) %>%
  group_by(celltype) %>%
  select(celltype, clone_id) %>%
  dplyr::count(celltype, clone_id) %>% mutate(clone_size = n) %>%
  select(-n)
countSize
## # A tibble: 180 × 3
## # Groups:   celltype [6]
##    celltype        clone_id    clone_size
##    <chr>           <chr>            <int>
##  1 HCL1            BCR_122_142       3006
##  2 HCL2            BCR_181_7         2564
##  3 Healthy B-Cells BCR_100_64           1
##  4 Healthy B-Cells BCR_101_82           1
##  5 Healthy B-Cells BCR_102_135          1
##  6 Healthy B-Cells BCR_103_129          1
##  7 Healthy B-Cells BCR_104_54           1
##  8 Healthy B-Cells BCR_105_46           1
##  9 Healthy B-Cells BCR_106_163          1
## 10 Healthy B-Cells BCR_107_34           1
## # ℹ 170 more rows
###################
# Donut charts of clonal composition
###################
## Calculate the size of each clone - don't subset to only cells with clone_id
# otherwise NA T-cells will not be in the plot later.
cloneCount <- plotTab.vdj %>% 
  group_by(celltype) %>%
  dplyr::count(clone_id) %>% ungroup()
cloneCount
## # A tibble: 186 × 3
##    celltype        clone_id        n
##    <chr>           <chr>       <int>
##  1 HCL1            BCR_122_142  3006
##  2 HCL1            <NA>          626
##  3 HCL2            BCR_181_7    2564
##  4 HCL2            <NA>          513
##  5 Healthy B-Cells BCR_100_64      1
##  6 Healthy B-Cells BCR_101_82      1
##  7 Healthy B-Cells BCR_102_135     1
##  8 Healthy B-Cells BCR_103_129     1
##  9 Healthy B-Cells BCR_104_54      1
## 10 Healthy B-Cells BCR_105_46      1
## # ℹ 176 more rows
cloneNum <- lapply(unique(plotTab.vdj$celltype), function(x) {
  subTab <- filter(plotTab.vdj, celltype == x) # subset to each patient
  countTab <- dplyr::count(subTab, clone_id)
  nclones <- countTab %>% nrow() # get the number of clones
  ncells_bcr <- nrow(subTab) # calculate the number of B or T cells per patient - important to subset plotTab to only B or T-cell clusters before doing this
  res <- data.frame(patID = x, nclones = nclones, ncells_bcr = ncells_bcr)
  res$type <- opt$subType
  res
}) %>% bind_rows()
cloneNum
##             patID nclones ncells_bcr
## 1            MCL1       2       3084
## 2            HCL1       2       3632
## 3            HCL2       2       3077
## 4            MCL3       2       3799
## 5            MCL2       2       2082
## 6 Healthy B-Cells     176        218
df <- left_join(cloneCount, cloneNum, by = c("celltype" = "patID")) %>%
  group_by(celltype) %>% mutate(clone_ratio = n / ncells_bcr) %>% # divide the size of the clone 
  select(celltype, clone_id, clone_ratio, nclones) %>% unique() %>%
  mutate(ymax = cumsum(clone_ratio), ymin = c(0, head(ymax, n= -1)))
df
## # A tibble: 186 × 6
## # Groups:   celltype [6]
##    celltype        clone_id    clone_ratio nclones    ymax    ymin
##    <chr>           <chr>             <dbl>   <int>   <dbl>   <dbl>
##  1 HCL1            BCR_122_142     0.828         2 0.828   0      
##  2 HCL1            <NA>            0.172         2 1       0.828  
##  3 HCL2            BCR_181_7       0.833         2 0.833   0      
##  4 HCL2            <NA>            0.167         2 1       0.833  
##  5 Healthy B-Cells BCR_100_64      0.00459     176 0.00459 0      
##  6 Healthy B-Cells BCR_101_82      0.00459     176 0.00917 0.00459
##  7 Healthy B-Cells BCR_102_135     0.00459     176 0.0138  0.00917
##  8 Healthy B-Cells BCR_103_129     0.00459     176 0.0183  0.0138 
##  9 Healthy B-Cells BCR_104_54      0.00459     176 0.0229  0.0183 
## 10 Healthy B-Cells BCR_105_46      0.00459     176 0.0275  0.0229 
## # ℹ 176 more rows
#### Overview over all patients/celltypes
pClone.Pat <- df %>% 
  ggplot(aes(ymax = ymax, ymin = ymin, xmax=4, xmin=3, fill = clone_id)) +
  geom_rect(colour = "black", size = 0.25) +
  coord_polar(theta="y") +
  xlim(c(2, 4)) +
  # ggtitle(paste0("Clonal Composition per Patient - ", opt$subType)) +
  facet_wrap(. ~ celltype, nrow =  1) +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5),
        strip.background = element_blank(),
        strip.text.x = element_text(size = 12.5),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA),
        legend.position = "none")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
pClone.Pat

#### Output individual patients
subPat <- c("HCL1", "HCL2", "MCL1", "MCL2", "MCL3", "Healthy B-Cells")

lapply(subPat, function(subPat) {
  pClone <- df %>% filter(celltype == subPat) %>%
  ggplot(aes(ymax = ymax, ymin = ymin, xmax=4, xmin=3, fill = clone_id)) +
  geom_rect(colour = "black", size = 0.25) +
  coord_polar(theta="y") +
  xlim(c(2, 4)) +
  ggtitle(paste0(subPat)) +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 75),
        legend.position = "none")
  pClone
  
  ggsave(plot = pClone, filename = paste0(opt$plot, "bcr_comp_", subPat, ".png"), 
         height = 20, width = 20, units = "cm")
})
## [[1]]
## [1] "plots/SFig4_5/bcr_comp_HCL1.png"
## 
## [[2]]
## [1] "plots/SFig4_5/bcr_comp_HCL2.png"
## 
## [[3]]
## [1] "plots/SFig4_5/bcr_comp_MCL1.png"
## 
## [[4]]
## [1] "plots/SFig4_5/bcr_comp_MCL2.png"
## 
## [[5]]
## [1] "plots/SFig4_5/bcr_comp_MCL3.png"
## 
## [[6]]
## [1] "plots/SFig4_5/bcr_comp_Healthy B-Cells.png"

Per patient UMAPs

patSub <- names(sce.pat)
sce.pat <- mclapply(patSub, mc.cores = ncores, function(x) {
  sce <- sce.pat[[x]]

  set.seed(100)
  sce <- runUMAP(sce, n_neighbors = 15, min_dist = 0.35, name = "UMAP", BPPARAM = mcparam)

})
names(sce.pat) <- patSub

patSub <- c("H371", "H431", "H279")

lapply(patSub, function(x) {
  sce <- sce.pat[[x]]
  
  set.seed(100)
  sce <- sce[,sample(1:ncol(sce), ncol(sce))]
  plotTab <- ggcells(sce)[[1]]
  
  if(x == "H371") {tlt <- "T-PLL1"
  } else if(x == "H431") {tlt <- "T-PLL2"
  } else if(x == "H279") {tlt <- "T-PLL3"} else {
    tlt <- x
  }  
  
  p <- plotTab %>% 
    ggplot(aes(x = UMAP.1, y = UMAP.2, fill = Treatment)) +
    geom_point(shape = 21, size = 3) +
    ggtitle(tlt) +
    lgd +
    theme_void() +
    theme(axis.text = element_blank(), 
          axis.ticks = element_blank(), 
          legend.position = "none", 
          plot.title = element_text(size = 20, hjust = 0.5)) +
    scale_fill_manual(values = c("DMSO" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000", 
                                 "Ibrutinib" = "#283593", "Everolimus" = "#43A047", 
                                 "Selumetinib" = "#F8BBD0"))
  
  ggsave(plot = p, filename = paste0(opt$plot, x, "_perpat_umap.png"), 
         height = 12.5, width = 12.5, units = "cm")
  
  p
})
## [[1]]

## 
## [[2]]

## 
## [[3]]

Output session info

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
##  [1] parallel  grid      stats4    stats     graphics  grDevices datasets 
##  [8] utils     methods   base     
## 
## other attached packages:
##  [1] sna_2.7-2                   statnet.common_4.9.0       
##  [3] scales_1.3.0                numbat_1.4.0               
##  [5] Matrix_1.6-1.1              airr_1.5.0                 
##  [7] ggsci_3.2.0                 BiocParallel_1.36.0        
##  [9] ggnet_0.1.0                 network_1.18.2             
## [11] tftargets_1.3               ggrepel_0.9.5              
## [13] readxl_1.4.3                edgeR_4.0.16               
## [15] limma_3.58.1                circlize_0.4.16            
## [17] ComplexHeatmap_2.18.0       gridExtra_2.3              
## [19] gprofiler2_0.2.3            scater_1.30.1              
## [21] scran_1.30.2                scuttle_1.12.0             
## [23] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
## [25] Biobase_2.62.0              GenomicRanges_1.54.1       
## [27] GenomeInfoDb_1.38.8         IRanges_2.36.0             
## [29] S4Vectors_0.40.2            BiocGenerics_0.48.1        
## [31] MatrixGenerics_1.14.0       matrixStats_1.3.0          
## [33] lubridate_1.9.3             forcats_1.0.0              
## [35] stringr_1.5.1               dplyr_1.1.4                
## [37] purrr_1.0.2                 readr_2.1.5                
## [39] tidyr_1.3.1                 tibble_3.2.1               
## [41] ggplot2_3.5.1               tidyverse_2.0.0            
## 
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22          splines_4.3.2            
##   [3] bitops_1.0-7              ggplotify_0.1.2          
##   [5] cellranger_1.1.0          polyclip_1.10-7          
##   [7] lifecycle_1.0.4           doParallel_1.0.17        
##   [9] lattice_0.21-9            MASS_7.3-60              
##  [11] magrittr_2.0.3            plotly_4.10.4            
##  [13] sass_0.4.9                rmarkdown_2.27           
##  [15] jquerylib_0.1.4           yaml_2.3.9               
##  [17] metapod_1.10.1            RColorBrewer_1.1-3       
##  [19] abind_1.4-5               zlibbioc_1.48.2          
##  [21] quadprog_1.5-8            hahmmr_1.0.0             
##  [23] ggraph_2.2.1              RCurl_1.98-1.14          
##  [25] yulab.utils_0.1.5         tweenr_2.0.3             
##  [27] GenomeInfoDbData_1.2.11   irlba_2.3.5.1            
##  [29] tidytree_0.4.6            dqrng_0.4.1              
##  [31] DelayedMatrixStats_1.24.0 codetools_0.2-19         
##  [33] DelayedArray_0.28.0       ggforce_0.4.2            
##  [35] tidyselect_1.2.1          shape_1.4.6.1            
##  [37] aplot_0.2.3               farver_2.1.2             
##  [39] ScaledMatrix_1.10.0       viridis_0.6.5            
##  [41] jsonlite_1.8.8            GetoptLong_1.0.5         
##  [43] BiocNeighbors_1.20.2      tidygraph_1.3.1          
##  [45] iterators_1.0.14          systemfonts_1.1.0        
##  [47] foreach_1.5.2             tools_4.3.2              
##  [49] ragg_1.3.2                treeio_1.29.0.002        
##  [51] Rcpp_1.0.12               glue_1.7.0               
##  [53] SparseArray_1.2.4         mgcv_1.9-0               
##  [55] xfun_0.45                 withr_3.0.0              
##  [57] fastmap_1.2.0             bluster_1.12.0           
##  [59] fansi_1.0.6               digest_0.6.36            
##  [61] rsvd_1.0.5                parallelDist_0.2.6       
##  [63] timechange_0.3.0          R6_2.5.1                 
##  [65] gridGraphics_0.5-1        textshaping_0.4.0        
##  [67] colorspace_2.1-0          Cairo_1.6-2              
##  [69] RhpcBLASctl_0.23-42       utf8_1.2.4               
##  [71] generics_0.1.3            renv_1.0.7               
##  [73] data.table_1.15.4         graphlayouts_1.1.1       
##  [75] httr_1.4.7                htmlwidgets_1.6.4        
##  [77] S4Arrays_1.2.1            uwot_0.2.2               
##  [79] pkgconfig_2.0.3           gtable_0.3.5             
##  [81] XVector_0.42.0            htmltools_0.5.8.1        
##  [83] clue_0.3-65               png_0.1-8                
##  [85] ggfun_0.1.5               knitr_1.48               
##  [87] rstudioapi_0.16.0         tzdb_0.4.0               
##  [89] rjson_0.2.21              coda_0.19-4.1            
##  [91] nlme_3.1-163              cachem_1.1.0             
##  [93] GlobalOptions_0.1.2       vipor_0.4.7              
##  [95] pillar_1.9.0              logger_0.3.0             
##  [97] vctrs_0.6.5               BiocSingular_1.18.0      
##  [99] beachmat_2.18.1           cluster_2.1.4            
## [101] beeswarm_0.4.0            evaluate_0.24.0          
## [103] cli_3.6.3                 locfit_1.5-9.10          
## [105] compiler_4.3.2            rlang_1.1.4              
## [107] crayon_1.5.3              labeling_0.4.3           
## [109] fs_1.6.4                  ggbeeswarm_0.7.2         
## [111] stringi_1.8.4             viridisLite_0.4.2        
## [113] munsell_0.5.1             lazyeval_0.2.2           
## [115] hms_1.1.3                 patchwork_1.2.0          
## [117] sparseMatrixStats_1.14.0  statmod_1.5.0            
## [119] highr_0.11                igraph_2.0.3             
## [121] memoise_2.0.1             scistreer_1.2.0          
## [123] RcppParallel_5.1.8        bslib_0.7.0              
## [125] phangorn_2.11.1           ggtree_3.10.1            
## [127] fastmatch_1.1-4           ape_5.8